#include "Jacobi.h"
#include "math.h"

T* diagonalized(CRSMatrix A)
{
    T* res=new T[A.numRow];
    for(unsigned int i=0; i<A.numRow; i++)
    {
        T tmp=0;
        unsigned int j=A.pointerB[i];
        while(A.column[j]<i)
        { j++; }
        if(A.column[j]==i)
            tmp=A.val[j];
        res[i]=tmp;
    }
    return res;
}

T deviation(T* x0, T* x, unsigned int n)
{
    T max=fabs(x0[0]-x[0]);
    for(unsigned int i=1; i<n; i++)
    {
        T tmp=fabs(x0[i]-x[i]);
        if(max<tmp) max=tmp;
    }
    return max;
}

unsigned int MethodJacobi(CRSMatrix A, T* x0, T* b, T* x, T E, T* dev)
{
    T* diag=diagonalized(A);
    T* Ax=Mul(A,x0);
    unsigned int iteration=1;
    for(unsigned int i=0; i<A.numRow; i++)
    {
        x[i]=x0[i]-(Ax[i]-b[i])/diag[i];
    }
    dev[0]=deviation(x0,x,A.numRow);
    while(dev[0]>E)
    {
        iteration++;
        delete(Ax);
        for(unsigned int i=0; i<A.numRow; i++)
        {
            x0[i]=x[i];
        }
        Ax=Mul(A, x);
        for(unsigned int i=0; i<A.numRow; i++)
        {
            x[i]=x[i]-(Ax[i]-b[i])/diag[i];
        }
        dev[0]=deviation(x0,x,A.numRow);
    }
    return iteration;
}